library(sf)
library(readxl)
library(dplyr)
library(plyr)
library(ggplot2)
library(afrihealthsites)
library(ggpubr)
library(afriadmin)
library(tmap)
library(cowplot)
library(forcats)
library(colorspace)
# Install Malawi MFL
malawi_MFL = read_excel("~/malawi-health-facilities-1/MHFR_Facilities 1.xlsx")
# Convert to sf
## omit NA's
new_malawi_MFL = na.omit(malawi_MFL)
## check for NA
any(is.na(new_malawi_MFL))
## [1] FALSE
## transform geometry columns into numeric
sapply(new_malawi_MFL, class)
## CODE NAME COMMON NAME OWNERSHIP TYPE STATUS
## "character" "character" "character" "character" "character" "character"
## ZONE DISTRICT DATE OPENED LATITUDE LONGITUDE
## "character" "character" "character" "character" "character"
new_malawi_MFL = transform(new_malawi_MFL, LATITUDE = as.numeric(LATITUDE),
LONGITUDE = as.numeric(LONGITUDE))
## Warning in eval(substitute(list(...)), `_data`, parent.frame()): NAs introduced
## by coercion
any(is.na(new_malawi_MFL)) ## check for NA
## [1] TRUE
new_malawi_MFL = na.omit(new_malawi_MFL) ## and omit
## convert to sf object
malawi_facilities_MFL = st_as_sf(new_malawi_MFL, coords = c("LONGITUDE", "LATITUDE"), dim = "XY")
malawi_facilities_MFL = st_set_crs(malawi_facilities_MFL, 4326) ## set CRS, is WGS84 right?
head(malawi_facilities_MFL)
## Simple feature collection with 6 features and 9 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 33.74129 ymin: -15.84 xmax: 35.09 ymax: -13.79742
## geographic CRS: WGS 84
## CODE NAME COMMON.NAME
## 1 MC010002 A + A private clinic A+A
## 2 BT240003 A-C Opticals A.C Opticals
## 3 BT240005 Akwezeke PVT Clinic Akwezeke Pvt
## 4 BT240006 AB Medical Clinic Abowa
## 5 LL040007 ABC Comm. Hospital ABC Clinic
## 6 LL040010 Achikondi Women Community Friendly Services Clinic Achikondi
## OWNERSHIP TYPE STATUS
## 1 Private Clinic Functional
## 2 Private Clinic Functional
## 3 Private Clinic Functional
## 4 Private Clinic Functional
## 5 Christian Health Association of Malawi (CHAM) Hospital Functional
## 6 Private Dispensary Functional
## ZONE DISTRICT DATE.OPENED geometry
## 1 Centrals West Zone Mchinji Jan 1st 75 POINT (33.88563 -13.79742)
## 2 South East Zone Blantyre Jan 1st 75 POINT (35.03 -15.8)
## 3 South East Zone Blantyre Jan 1st 75 POINT (35.09 -15.84)
## 4 South East Zone Blantyre Jan 1st 75 POINT (35.09 -15.84)
## 5 Centrals West Zone Lilongwe Jan 1st 75 POINT (33.74129 -13.96816)
## 6 Centrals West Zone Lilongwe Jan 1st 75 POINT (33.7793 -13.95473)
Overview:
# Re-order the facility types
malawi_MFL$TYPE = as.factor(malawi_MFL$TYPE)
malawi_MFL$TYPE = factor(malawi_MFL$TYPE, levels = c("Central Hospital", "District Hospital", "Hospital", "Health Centre", "Clinic", "Health Post", "Dispensary", "Private", "Unclassified"))
# Number of each type + ownership
facility_types_MFL = as.data.frame(table(malawi_MFL$TYPE, malawi_MFL$OWNERSHIP))
head(facility_types_MFL)
## Var1 Var2 Freq
## 1 Central Hospital Aquaid Lifeline 0
## 2 District Hospital Aquaid Lifeline 0
## 3 Hospital Aquaid Lifeline 0
## 4 Health Centre Aquaid Lifeline 0
## 5 Clinic Aquaid Lifeline 0
## 6 Health Post Aquaid Lifeline 0
levels(facility_types_MFL$Var2)
## [1] "Aquaid Lifeline"
## [2] "Christian Health Association of Malawi (CHAM)"
## [3] "Government"
## [4] "Mission/Faith-based (other than CHAM)"
## [5] "Non-Government"
## [6] "Other"
## [7] "Parastatal"
## [8] "Private"
# Only showing government, private and CHAM ownership
facility_types_MFL$new_Var2 = revalue(facility_types_MFL$Var2, c("Aquaid Lifeline"="Other", "Non-Government"="Other", "Parastatal"="Other", "Mission/Faith-based (other than CHAM)"="Other"))
facility_types_MFL$new_Var2 = factor(facility_types_MFL$new_Var2, levels = c("Government", "Private", "Christian Health Association of Malawi (CHAM)", "Other"))
## bar plot of no. of facility types
plot_facility_types_MFL = ggplot(facility_types_MFL, aes(x=Var1, y=Freq, fill=new_Var2)) + geom_bar(position = "stack", stat = "identity")
plot_facility_types_MFL = plot_facility_types_MFL + labs(x = "Facility types", y = "Frequency", fill="Ownership") + scale_fill_brewer(palette = "Set2") + coord_flip() + theme_minimal() + ggtitle("MFL") + theme(axis.title.y = element_blank(), legend.title = element_text(size = 8, face = "bold"), legend.text = element_text(size = 8), plot.title = element_text(face = "bold", hjust = 0), legend.position = "bottom")
plot_facility_types_MFL
par(mar=c(11,4,4,4))
# Re-order ownership
malawi_MFL$OWNERSHIP = as.factor(malawi_MFL$OWNERSHIP)
malawi_MFL$OWNERSHIP = factor(malawi_MFL$OWNERSHIP, levels = c("Government", "Private", "Christian Health Association of Malawi (CHAM)", "Non-Government", "Mission/Faith-based (other than CHAM)", "Other", "Parastatal", "Aquaid Lifeline"))
# Number of each type of ownership
ownership_MFL = as.data.frame(table(malawi_MFL$OWNERSHIP))
## bar plot of ownership
plot_ownership_MFL = ggplot(ownership_MFL, aes(x=forcats::fct_relabel(Var1,stringr::str_wrap,width = 16), y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_ownership_MFL = plot_ownership_MFL + labs(x="Ownership", y = "Frequency") + coord_flip() + theme_minimal() + ggtitle("MFL") + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold"))
plot_ownership_MFL
Focuses on facilities run by government, faith-based organisations, NGO’s and local authorities. Covers 50 countries in sub-Saharan Africa. Sources of information include health sector reports, websites run by national or international organisations and personal communications
If MFL was available it was used. More than one datasource was often used to compile facility list
Private facilities are excluded, duplicates removed, name errors corrected and name variations were matched. Missing info was added with the use of other datasources.
Now hosted by the WHO Global Malaria Programme, last update February 2019
Malawi datasources includes MFL, https://data.humdata.org/dataset/malawi-health and http://www.cham.org.mw/uploads/7/3/0/8/73088105/cham_health_facilities_-_1_june_2016.pdf
At time of publishing, 639 facilities with 9 missing coordinates, not been updated since
Data includes facility name, type, ownership, source of location and reclassified facility types
# Malawi WHO data.frame
malawi_WHO <- afrihealthsites("malawi", datasource='who', plot=FALSE, returnclass='dataframe')
head(malawi_WHO)
## # A tibble: 6 x 10
## Country Admin1 `Facility name` `Facility type` Ownership Lat Long
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 Malawi Centr… 80 Block Clinic Clinic MoH -12.9 33.4
## 2 Malawi Centr… ABC Community … Clinic FBO -14.0 33.7
## 3 Malawi Centr… Adventist Heal… Health Centre FBO -14.0 33.8
## 4 Malawi Centr… Alinafe Commun… Community Hosp… FBO -13.4 34.2
## 5 Malawi Centr… Area 18 Health… Health Centre MoH -13.9 33.8
## 6 Malawi Centr… Area 25 Health… Health Centre MoH -13.9 33.8
## # … with 3 more variables: `LL source` <chr>, iso3c <chr>,
## # facility_type_9 <chr>
# Re order facility types and ownership
malawi_WHO$`Facility type` = as.factor(malawi_WHO$`Facility type`)
malawi_WHO$`Facility type` = factor(malawi_WHO$`Facility type`, levels = c("Central Hospital", "District Hospital", "Mission Hospital", "Rural Hospital", "Community Hospital", "Health Centre", "Clinic", "Health Post/Dispensary"))
malawi_WHO$Ownership = as.factor(malawi_WHO$Ownership)
malawi_WHO$Ownership = factor(malawi_WHO$Ownership, levels = c("MoH", "FBO", "Local authority", "NGO"))
# No. of original facility types + ownership
facility_types_WHO = as.data.frame(table(malawi_WHO$`Facility type`, malawi_WHO$Ownership))
head(facility_types_WHO)
## Var1 Var2 Freq
## 1 Central Hospital MoH 4
## 2 District Hospital MoH 24
## 3 Mission Hospital MoH 0
## 4 Rural Hospital MoH 17
## 5 Community Hospital MoH 0
## 6 Health Centre MoH 332
## bar plot of original facility types
plot_facility_types_WHO = ggplot(facility_types_WHO, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(position = "stack", stat = "identity")
plot_facility_types_WHO = plot_facility_types_WHO + labs(x = "Facility types", y = "Frequency", fill="Ownership") + scale_fill_brewer(palette = "Set2") + coord_flip() + theme_minimal() + ggtitle("WHO") + theme(axis.title.y = element_blank(), legend.title = element_text(size = 8, face = "bold"), legend.text = element_text(size = 8), plot.title = element_text(face = "bold", hjust = 0), legend.position = "bottom") + expand_limits(y=c(0,600))
plot_facility_types_WHO
# Re order
malawi_WHO$facility_type_9 = as.factor(malawi_WHO$facility_type_9)
malawi_WHO$facility_type_9 = factor(malawi_WHO$facility_type_9, levels = c("Hospital", "Health Centre", "Health Clinic", "Health Post", "Community Health Unit"))
# No. of reclassified facility types
RC_facility_types_WHO = as.data.frame(table(malawi_WHO$facility_type_9))
RC_facility_types_WHO
## Var1 Freq
## 1 Hospital 80
## 2 Health Centre 457
## 3 Health Clinic 22
## 4 Health Post 87
## 5 Community Health Unit 2
## bar plot of reclassified facility types
plot_RC_facility_types_WHO = ggplot(RC_facility_types_WHO, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_RC_facility_types_WHO = plot_RC_facility_types_WHO + labs(x = "Reclassified facility types", y = "Frequency") + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold")) + ggtitle("WHO") + expand_limits(y=c(0,600))
plot_RC_facility_types_WHO
# Types of ownership
ownership_WHO = as.data.frame(table(malawi_WHO$Ownership))
ownership_WHO
## Var1 Freq
## 1 MoH 467
## 2 FBO 173
## 3 Local authority 5
## 4 NGO 3
## bar plot of ownership
plot_ownership_WHO = ggplot(ownership_WHO, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_ownership_WHO = plot_ownership_WHO + labs(x="Ownership", y = "Frequency") + coord_flip() + theme_minimal() + ggtitle("WHO") + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold")) + expand_limits(y=c(0,600))
plot_ownership_WHO
Both data sources contain no information on services available, capacity or equipment. MFL does state whether facility is functional.
Classification of MFL facilities aligns more with the structure of the health care system in Malawi (community, primary, secondary, tertiary), it differentiates central hospitals from district and other hospitals. WHO has additional rural and mission hospitals, where do they fit in?
https://www.health.gov.mw/index.php/2016-01-06-19-58-23/national-aids states that at community level, health posts, dispensaries and maternity clinics offer services. Primary includes health centers and community hospitals, secondary consists of district and some CHAM hospitals, tertiary includes central hospitals.
Analysis:
## tmap mode set to interactive viewing
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
Qs to address?:
# Qs 1 - how many intersect?
## convert malawi_WHO to sf object
class(malawi_WHO)
## [1] "tbl_df" "tbl" "data.frame"
any(is.na(malawi_WHO))
## [1] TRUE
new_malawi_WHO = na.omit(malawi_WHO) ## omit NA
sf_malawi_WHO = st_as_sf(new_malawi_WHO, coords = c("Long", "Lat"), dim = "XY")
sf_malawi_WHO = st_set_crs(sf_malawi_WHO, 4326)
## st_intersection
intersect_WHO_MFL = st_intersection(x=sf_malawi_WHO, y=malawi_facilities_MFL)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
intersect_WHO_MFL ## only 2 intersect directly, so are same up to 5 decimal places?
## Simple feature collection with 2 features and 17 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 33.29456 ymin: -11.53894 xmax: 33.41925 ymax: -11.45836
## geographic CRS: WGS 84
## # A tibble: 2 x 18
## Country Admin1 Facility.name Facility.type Ownership LL.source iso3c
## * <chr> <chr> <chr> <fct> <fct> <chr> <chr>
## 1 Malawi North… Euthini Heal… Health Centre MoH GPS MWI
## 2 Malawi North… Madede Healt… Health Centre MoH GPS MWI
## # … with 11 more variables: facility_type_9 <fct>, CODE <chr>, NAME <chr>,
## # COMMON.NAME <chr>, OWNERSHIP <chr>, TYPE <fct>, STATUS <chr>, ZONE <chr>,
## # DISTRICT <chr>, DATE.OPENED <chr>, geometry <POINT [°]>
# Try merge_points()
merge_points("malawi", datasources=c('who', 'healthsites'), dist_same_m=50)
## although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
## malawi num points: who : 639 healthsites : 266 shared at dist thresh 50 m : 52 after merging: 853
# Facilities per region
# MFL, admin1
facility_admin1_MFL = st_intersects(malawi_admin1, malawi_facilities_MFL, sparse = TRUE)
malawi_admin1$facility = lengths(facility_admin1_MFL)
head(malawi_admin1) # gave no. of facilities per region
## Simple feature collection with 3 features and 6 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 32.67036 ymin: -17.12952 xmax: 35.91857 ymax: -9.368326
## geographic CRS: WGS 84
## shapeName shapeISO shapeID shapeGroup shapeType
## 1 Central Region MW-C MWI-ADM1-3_0_0-B1 MWI ADM1
## 2 Southern Region MW-S MWI-ADM1-3_0_0-B2 MWI ADM1
## 3 Northern Region MW-N MWI-ADM1-3_0_0-B3 MWI ADM1
## geometry facility
## 1 MULTIPOLYGON (((33.5486 -12... 508
## 2 MULTIPOLYGON (((34.58346 -1... 609
## 3 MULTIPOLYGON (((34.02092 -1... 248
# Intersect opposite to get admin1 region for each facility
admin1_facility_MFL = unlist(st_intersects(malawi_facilities_MFL, malawi_admin1, sparse = TRUE))
# error when adding admin1_facility_MFL data to malawi_facilities_MFL,
# no. of rows are less in admin1_, why?
# trying to figure out which points don't intersect
intersect_admin1_MFL = st_intersection(malawi_facilities_MFL, malawi_admin1) # returns data frame of matches
intersect_admin1_MFL = intersect_admin1_MFL[-c(10:16)] # remove columns
# adding region data for admin1 to intersect_admin1_MFL (using admin1_facility_MFL from previously)
intersect_admin1_MFL$region_admin1 = admin1_facility_MFL
# working out which facilities did not intersect
intersect_admin1_MFL = as.data.frame(intersect_admin1_MFL) # convert to data frame to use in anti_join
head(intersect_admin1_MFL)
## CODE NAME
## 1 MC010002 A + A private clinic
## 5 LL040007 ABC Comm. Hospital
## 6 LL040010 Achikondi Women Community Friendly Services Clinic
## 7 LL040011 Lilongwe Adventist Hospital
## 8 LL040012 Adventist Health Centre Area 15
## 9 LL040013 Africa Leaf Clinic Kanengo
## COMMON.NAME OWNERSHIP
## 1 A+A Private
## 5 ABC Clinic Christian Health Association of Malawi (CHAM)
## 6 Achikondi Private
## 7 Adventist Health Centre Christian Health Association of Malawi (CHAM)
## 8 Adventist Christian Health Association of Malawi (CHAM)
## 9 Africa Leaf Clinic Kanengo Private
## TYPE STATUS ZONE DISTRICT DATE.OPENED
## 1 Clinic Functional Centrals West Zone Mchinji Jan 1st 75
## 5 Hospital Functional Centrals West Zone Lilongwe Jan 1st 75
## 6 Dispensary Functional Centrals West Zone Lilongwe Jan 1st 75
## 7 Hospital Functional Centrals West Zone Lilongwe Jan 1st 83
## 8 Health Centre Functional Centrals West Zone Lilongwe Jan 1st 75
## 9 Clinic Non-functional Centrals West Zone Lilongwe Jan 1st 75
## geometry region_admin1
## 1 POINT (33.88563 -13.79742) 1
## 5 POINT (33.74129 -13.96816) 2
## 6 POINT (33.7793 -13.95473) 2
## 7 POINT (33.7793 -13.95473) 2
## 8 POINT (33.7793 -13.95473) 1
## 9 POINT (33.80487 -13.8898) 1
no_intersect_admin1_MFL = anti_join(malawi_facilities_MFL, intersect_admin1_MFL) # shows which did not intersect
# input 0s in the region_admin1 column
no_intersect_admin1_MFL$region_admin1 = rep.int(0, times = 61)
# turn back into sf
no_intersect_admin1_MFL = st_as_sf(no_intersect_admin1_MFL)
intersect_admin1_MFL = st_as_sf(intersect_admin1_MFL)
# bind data frames together to get new data frame with admin1 regions assigned correctly
malawi_facilities_MFL_admin1 = rbind(intersect_admin1_MFL, no_intersect_admin1_MFL)
tail(malawi_facilities_MFL_admin1)
## Simple feature collection with 6 features and 10 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 1 ymin: -33.45 xmax: 26.56 ymax: -1
## geographic CRS: WGS 84
## CODE NAME COMMON.NAME
## 5610 MN281537 GOOD SAMARITAN PVT CLINIC GOOD SAMARITAN PVT CLINIC
## 5710 MC011538 Kazyzyo Dispensary Kazyozyo Dispensary
## 5810 MH201539 Cape Maclear Chembe
## 5910 MH201540 Nankhwali Nankhwali
## 6010 MH201541 Kapire Health Centre Kapire
## 6110 ZA231548 Prime PVT CLINIC Jali pvt
## OWNERSHIP TYPE
## 5610 Private Private
## 5710 Government Dispensary
## 5810 Government Health Post
## 5910 Christian Health Association of Malawi (CHAM) Health Centre
## 6010 Christian Health Association of Malawi (CHAM) Health Centre
## 6110 Private Private
## STATUS ZONE DISTRICT
## 5610 Functional South East Zone Mwanza
## 5710 Functional Centrals West Zone Mchinji
## 5810 Pending Operation (Under construction) South West Zone Mangochi
## 5910 Functional South West Zone Mangochi
## 6010 Functional South West Zone Mangochi
## 6110 Functional South West Zone Zomba
## DATE.OPENED region_admin1 geometry
## 5610 Jan 1st 75 0 POINT (3 -2)
## 5710 Jan 1st 90 0 POINT (26.56 -33.45)
## 5810 Jan 1st 75 0 POINT (2 -1)
## 5910 Jan 1st 75 0 POINT (2 -1)
## 6010 Jan 1st 75 0 POINT (1 -1)
## 6110 Jan 1st 75 0 POINT (1 -1)
# !!! something is wrong, facilities with same district show different
# value for region. Was data for district inputted incorrectly or are
# coordinates wrong? Carry on for now.
# use new data frame for plots
# Map with admin1 layer
tmap_admin1_MFL = tmap_facilities_MFL + tm_shape(st_geometry(malawi_admin1)) + tm_borders()
tmap_admin1_MFL
## Map - number of facilities per region
tmap_admin1_MFL2 = tm_shape(st_geometry(malawi_admin1)) + tm_borders() + tm_shape(malawi_admin1) + tm_fill("facility", style = "cat", palette = sequential_hcl(4, palette = "YlGnBu", rev = TRUE), title = "MFL facilities")
tmap_admin1_MFL2
# Plot of facilities in each admin1 region
## transform columns
malawi_facilities_MFL_admin1$region_admin1 = as.factor(malawi_facilities_MFL_admin1$region_admin1)
malawi_facilities_MFL_admin1$TYPE = as.factor(malawi_facilities_MFL_admin1$TYPE)
## freq of facility types by region
df_facility_admin1_MFL = as.data.frame(table(malawi_facilities_MFL_admin1$region_admin1,
malawi_facilities_MFL_admin1$TYPE))
head(df_facility_admin1_MFL)
## Var1 Var2 Freq
## 1 0 Central Hospital 0
## 2 1 Central Hospital 0
## 3 2 Central Hospital 4
## 4 3 Central Hospital 0
## 5 0 District Hospital 1
## 6 1 District Hospital 7
## subsetting by region for plots
# region 1 = Central
admin1_region1 = filter(df_facility_admin1_MFL, Var1 == 1)
admin1_region1_plot = ggplot(admin1_region1, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "MFL Central Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 250))
# region 2 = Southern
admin1_region2 = filter(df_facility_admin1_MFL, Var1 == 2)
admin1_region2_plot = ggplot(admin1_region2, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "MFL Southern Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0,250))
# region 3 = Northern
admin1_region3 = filter(df_facility_admin1_MFL, Var1 == 3)
admin1_region3_plot = ggplot(admin1_region3, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "MFL Northern Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 150))
# MFL, admin2
facility_admin2_MFL = st_intersects(malawi_admin2, malawi_facilities_MFL, sparse = TRUE)
malawi_admin2$facility = lengths(facility_admin2_MFL)
sum(malawi_admin2$facility) # also has 1365 intersections
## [1] 1365
# Map - number of facilities per district
tmap_admin2_MFL = tm_shape(st_geometry(malawi_admin2)) + tm_borders() + tm_shape(malawi_admin2) + tm_fill("facility", palette = sequential_hcl(6, palette = "YlGnBu", rev = TRUE), title = "MFL facilities") + tm_layout(title = "MFL")
tmap_admin2_MFL
# WHO, admin1
# Map with admin1 layer
tmap_admin1_WHO = tmap_facilities_WHO + tm_shape(st_geometry(malawi_admin1)) + tm_borders()
tmap_admin1_WHO
## intersecting for Map
facility_admin1_WHO = st_intersects(malawi_admin1, sf_malawi_WHO, sparse = TRUE)
malawi_admin1$facility_WHO = lengths(facility_admin1_WHO)
head(malawi_admin1)
## Simple feature collection with 3 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 32.67036 ymin: -17.12952 xmax: 35.91857 ymax: -9.368326
## geographic CRS: WGS 84
## shapeName shapeISO shapeID shapeGroup shapeType
## 1 Central Region MW-C MWI-ADM1-3_0_0-B1 MWI ADM1
## 2 Southern Region MW-S MWI-ADM1-3_0_0-B2 MWI ADM1
## 3 Northern Region MW-N MWI-ADM1-3_0_0-B3 MWI ADM1
## geometry facility facility_WHO
## 1 MULTIPOLYGON (((33.5486 -12... 508 238
## 2 MULTIPOLYGON (((34.58346 -1... 609 271
## 3 MULTIPOLYGON (((34.02092 -1... 248 130
## Map - number of facilities per region
tmap_admin1_WHO2 = tm_shape(st_geometry(malawi_admin1)) + tm_borders() + tm_shape(malawi_admin1) + tm_fill("facility_WHO", style = "cat", title = "WHO facilities", palette = sequential_hcl(4, palette = "YlGnBu", rev = TRUE))
tmap_admin1_WHO2
# freq of facility types by region
malawi_WHO$Admin1 = as.factor(malawi_WHO$Admin1)
df_facility_admin1_WHO = as.data.frame(table(malawi_WHO$Admin1,
malawi_WHO$`Facility type`))
head(df_facility_admin1_WHO)
## Var1 Var2 Freq
## 1 Central Central Hospital 1
## 2 Northern Central Hospital 1
## 3 Southern Central Hospital 2
## 4 Central District Hospital 9
## 5 Northern District Hospital 6
## 6 Southern District Hospital 9
# Central
admin1_central = filter(df_facility_admin1_WHO, Var1 == "Central")
admin1_central_plot = ggplot(admin1_central, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "WHO Central Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 250))
# Southern
admin1_southern = filter(df_facility_admin1_WHO, Var1 == "Southern")
admin1_southern_plot = ggplot(admin1_southern, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "WHO Southern Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0,250))
# Northern
admin1_northern = filter(df_facility_admin1_WHO, Var1 == "Northern")
admin1_northern_plot = ggplot(admin1_northern, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "WHO Northern Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 150))
# WHO, admin2
facility_admin2_WHO = st_intersects(malawi_admin2, sf_malawi_WHO, sparse = TRUE)
malawi_admin2$facility_WHO = lengths(facility_admin2_WHO)
# Map - number of facilities per district
tmap_admin2_WHO = tm_shape(st_geometry(malawi_admin2)) + tm_borders() + tm_shape(malawi_admin2) + tm_fill("facility_WHO", palette = sequential_hcl(4, palette = "YlGnBu", rev = TRUE), title = "WHO facilities")
tmap_admin2_WHO
# final plots
plot_admin1_central = align_plots(admin1_region1_plot, admin1_central_plot, align = "v")
ggdraw(plot_admin1_central[[1]])
ggdraw(plot_admin1_central[[2]])
plot_admin1_southern = align_plots(admin1_region2_plot, admin1_southern_plot, align = "v")
ggdraw(plot_admin1_southern[[1]])
ggdraw(plot_admin1_southern[[2]])
plot_admin1_northern = align_plots(admin1_region3_plot, admin1_northern_plot, align = "v")
ggdraw(plot_admin1_northern[[1]])
ggdraw(plot_admin1_northern[[2]])
# Function that recognises different facility types
sort_facilities = function(data, facility_column) {
for(i in levels(facility_column)) {
if(grepl("Central", facility_column, fixed = TRUE) == TRUE){
output = grep("Central", facility_column, fixed = TRUE)
t4 = slice(data, c(output[]))
}
return(t4)
}
}